﻿using LxBasic;
using LxCore;

namespace LxCore
{
    public class SIMPLEVelocity : ZoneHead
    {
        //local
        public NSEquation U;
        public NSEquation V;
        public NSEquation W;
        protected bool IsVelInit { get => U.IsXInit; set => U.IsXInit = V.IsXInit = W.IsXInit = value; }
        protected bool IsGradInit { get => U.IsGradInit; set => U.IsGradInit = V.IsGradInit = W.IsGradInit = value; }

        //share
        public IField<double> UVW_Gamma;
        public IField<Vector3d> UVW_Source_NoPressure_b;
        public IField<Vector3d> UVW_Source_NoPressure_a;
        public IField<Vector3d> UVW_IniValue;
        public IField<double> DynamicViscosity;

        public SIMPLEFields InSimpleFields { get; set; }
        public VarArrBuilder3.DirectS GetInnCellVelVars() { return new VarArrBuilder3.DirectS(U.InnX, V.InnX, W.InnX); }
        public VarArrBuilder3.DirectS GetBouCellVelVars() { return new VarArrBuilder3.DirectS(U.BouX, V.BouX, W.BouX); }
        //public MatrixR3C3 GetVelGrad(int cellIndex)//与红皮书不一致，红皮书疑似有误
        //{
        //    Vector3d gradU = U_Grad[cellIndex];
        //    Vector3d gradV = V_Grad[cellIndex];
        //    Vector3d gradW = W_Grad[cellIndex];
        //    return new MatrixR3C3(
        //        gradU.X, gradU.Y, gradU.Z,
        //        gradV.X, gradV.Y, gradV.Z,
        //        gradW.X, gradW.Y, gradW.Z);
        //}
        //public MatrixR3C3 GetVelGradT(int cellIndex)
        //{
        //    Vector3d u = U_Grad[cellIndex];
        //    Vector3d v = V_Grad[cellIndex];
        //    Vector3d w = W_Grad[cellIndex];
        //    return new MatrixR3C3(
        //        u.X, v.X, w.X,
        //        u.Y, v.Y, w.Y,
        //        u.Z, v.Z, w.Z);
        //}

        public override void CheckCreate(FvZone zone)
        {
            base.CheckCreate(zone);

            CheckCreator.Instance(ref U);
            U.CheckCreate(zone);

            CheckCreator.Instance(ref V);
            V.CheckCreate(zone);

            CheckCreator.Instance(ref W);
            W.CheckCreate(zone);
        }
        public override void Reuse(ZoneHead oldOne)
        {
            base.Reuse(oldOne);
            if (oldOne is SIMPLEVelocity o)
            {
                this.U = o.U;
                this.V = o.V;
                this.W = o.W;
            }
        }
        #region 设置Vector
        public void SetZone_CommonTerm(IField<double> UVW_gamma, IField<double> dynamicViscosity, IField<Vector3d> UVW_source_noPressure_b, IField<Vector3d> UVW_source_noPressure_a, IField<Vector3d> UVW_iniValue, IParameter<double> relaxFactor)
        {
            this.UVW_Gamma = UVW_gamma;
            this.DynamicViscosity = dynamicViscosity;
            this.UVW_Source_NoPressure_b = UVW_source_noPressure_b;
            this.UVW_Source_NoPressure_a = UVW_source_noPressure_a;
            this.UVW_IniValue = UVW_iniValue;

            setUVWCommonTerm(U);
            setUVWCommonTerm(V);
            setUVWCommonTerm(W);
            void setUVWCommonTerm(NSEquation uvw)
            {
                uvw.SetZone_CommonTerm(ConstField.Double0, ConstField.Double0, null, relaxFactor);
                uvw.SetZone_DiffTerm(this.UVW_Gamma);
            }
        }
        public void SetZone_ConvTerm(ConvectionSchemes convectionSchemes)
        {
            U.SetZone_ConvTerm(ConstField.Double1, convectionSchemes, InSimpleFields.MassFlux);
            V.SetZone_ConvTerm(ConstField.Double1, convectionSchemes, InSimpleFields.MassFlux);
            W.SetZone_ConvTerm(ConstField.Double1, convectionSchemes, InSimpleFields.MassFlux);
        }

        #endregion
        #region 初始化相关
        public override void Initialize()
        {
            if (!IsVelInit)
            {
                MultiThreads.RunParallelTask((threadInfo) =>
                {
                    setVelDirect(UVW_IniValue, Zone.Inner, threadInfo);
                });//slipwall边界的初始化依赖内部速度场，因此分两步
                MultiThreads.RunParallelTask((threadInfo) =>
                {
                    var mbs = InSimpleFields.SimpleBounds;
                    for (int i1 = 0; i1 < mbs.Length; i1++)
                    {
                        switch (mbs[i1])
                        {
                            case SIMPLEFields.NoSlipWall noslip:
                                setVelOfNoSlipWall(noslip.BoundVelocity, Zone.Bounds[i1], threadInfo);
                                break;
                            case SIMPLEFields.VelocityInlet velInlet:
                                setVelDirect(velInlet.Velocity, Zone.Bounds[i1], threadInfo);
                                break;
                            case SIMPLEFields.Symmetry sym:
                                setVelOfSlipWall(Zone.Bounds[i1], threadInfo);
                                break;
                            case SIMPLEFields.PressureInOut p:
                                setVelDirect(UVW_IniValue, Zone.Bounds[i1], threadInfo);
                                break;
                        }
                    }
                });
                IsVelInit = true;
            }
            if (!IsGradInit)
            {
                MultiThreads.RunParallelTask((Action<ThreadInfo>)U.CalGradient_LSQ + V.CalGradient_LSQ + W.CalGradient_LSQ);
                IsGradInit = true;
            }
            U.TransData?.Initialize();
            V.TransData?.Initialize();
            W.TransData?.Initialize();
        }
        protected void setVelDirect(IField<Vector3d> UVWValue, IIndexHost indexHost, ThreadInfo threadInfo)
        {
            if (indexHost is FvInner)
            {
                set(U.InnX, V.InnX, W.InnX);
            }
            else
            {
                set(U.BouX, V.BouX, W.BouX);
            }
            void set(double[] X, double[] Y, double[] Z)
            {
                (int i1, int endp1) = indexHost.CellIdx.Split(threadInfo);
                if (UVWValue.IsParameter(out Vector3d cv))
                {
                    if (cv == Vector3d.Zero)
                    {
                        Array.Clear(X, i1, endp1 - i1);
                        Array.Clear(Y, i1, endp1 - i1);
                        Array.Clear(Z, i1, endp1 - i1);
                    }
                    else
                    {
                        Array.Fill(X, cv.X, i1, endp1 - i1);
                        Array.Fill(Y, cv.Y, i1, endp1 - i1);
                        Array.Fill(Z, cv.Z, i1, endp1 - i1);
                    }
                }
                else
                {
                    var vel = UVWValue.GetCellVarArr(indexHost);
                    for (; i1 < endp1; i1++)
                    {
                        Vector3d v = vel.GetCellVar(i1);
                        X[i1] = v.X;
                        Y[i1] = v.Y;
                        Z[i1] = v.Z;
                    }
                }
            }
        }
        protected void setVelOfNoSlipWall(IField<Vector3d> UVWValue, FvBound bound, ThreadInfo threadInfo)
        {
            if (UVWValue.IsEqual0())
            {
                setVelDirect(UVWValue, bound, threadInfo);
                return;
            }
            (int i1, int endp1) = bound.FaceIdx.Split(threadInfo);
            var vel = UVWValue.GetCellVarArr(bound);
            double[] X = U.BouX; double[] Y = V.BouX; double[] Z = W.BouX;
            for (; i1 < endp1; i1++)
            {
                ref var face = ref BouFaces[i1];
                Vector3d v = vel.GetCellVar(i1);
                v = Vector3d.ProjectVectorToPlane(v, face.Normal);
                X[i1] = v.X;
                Y[i1] = v.Y;
                Z[i1] = v.Z;
            }
        }
        protected void setVelOfSlipWall(FvBound bound, ThreadInfo threadInfo)
        {
            (int i1, int endp1) = bound.FaceIdx.Split(threadInfo);
            var innVel = GetInnCellVelVars();
            double[] X = U.BouX; double[] Y = V.BouX; double[] Z = W.BouX;
            for (; i1 < endp1; i1++)
            {
                ref var face = ref BouFaces[i1];
                Vector3d v = innVel.GetCellVar(face.OSideCell);
                v = Vector3d.ProjectVectorToPlane(v, face.Normal);
                X[i1] = v.X;
                Y[i1] = v.Y;
                Z[i1] = v.Z;
            }
        }
        protected void formBoundFactorOfPressreInOut(FvBound bound, SIMPLEFields.PressureInOut pbound, ThreadInfo threadInfo)
        {
            double[] bouU = U.BouX; double[] bouV = V.BouX; double[] bouW = W.BouX;
            double[] innU = U.InnX; double[] innV = V.InnX; double[] innW = W.InnX;
            Vector3d[] gradU = U.Grad; Vector3d[] gradV = V.Grad; Vector3d[] gradW = W.Grad;
            var factors = pbound.Factors;
            int start = bound.FaceIdx.Start;
            (int i1, int endp1) = bound.FaceIdx.Split(threadInfo);
            for (; i1 < endp1; i1++)
            {
                ref var f = ref BouFaces[i1];
                factors[i1 - start] = new Vector3d(
                    -f.CalNonOrthogonalityCorrectionOfBound(gradU, innU, bouU[i1]) / f.gDiff,
                    -f.CalNonOrthogonalityCorrectionOfBound(gradV, innV, bouV[i1]) / f.gDiff,
                    -f.CalNonOrthogonalityCorrectionOfBound(gradW, innW, bouW[i1]) / f.gDiff
                    );
            }
        }
        protected void updateBoundValueOfPressureInOut(FvBound bound, SIMPLEFields.PressureInOut pbound, ThreadInfo threadInfo)
        {
            (int i1, int endp1) = bound.FaceIdx.Split(threadInfo);
            var indir = pbound.Direction.GetFaceVarArr(bound);
            var pressure = pbound.Pressure.GetFaceVarArr(bound);
            var density = InSimpleFields.Density.GetFaceVarArr(bound);
            var innVels = GetInnCellVelVars();
            var InnDc = InSimpleFields.Pressure.InnDc;
            var innP = InSimpleFields.Pressure.InnX;
            var bouP = InSimpleFields.Pressure.BouX;
            var pGrad = InSimpleFields.Pressure.Grad;
            double[] bouU = U.BouX; double[] bouV = V.BouX; double[] bouW = W.BouX;
            double[] BouMassFlux = InSimpleFields.MassFlux.BouMassFlux;
            double[] MassImbalance = InSimpleFields.MassFlux.MassImbalance;
            var factors = pbound.Factors;
            int start = bound.FaceIdx.Start;
            for (; i1 < endp1; i1++)
            {
                ref var face = ref BouFaces[i1];
                //要区分流动倾向性
                double massflux = BouMassFlux[i1];
                double massImb = MassImbalance[face.OSideCell];
                Vector3d cV = innVels.GetCellVar(face.OSideCell);
                if (massflux > 0)////有明确的流出时，确定为流出边界massflux - massImb > 0
                {
                    //
                    Vector3d bVOut = cV + factors[i1 - start];
                    if (Vector3d.DotProduct(bVOut, face.Normal) >= 0)
                    {
                        bouU[i1] = bVOut.X;
                        bouV[i1] = bVOut.Y;
                        bouW[i1] = bVOut.Z;
                        continue;
                    }
                    //
                    Vector3d dir = indir.GetFaceVar(i1);
                    Vector3d bVIn = Vector3d.DotProduct(cV, dir) / dir.LengthSquared() * dir;
                    if (Vector3d.DotProduct(bVIn, face.Normal) <= 0)
                    {
                        bouU[i1] = bVIn.X;
                        bouV[i1] = bVIn.Y;
                        bouW[i1] = bVIn.Z;
                        continue;
                    }
                }
                else
                {
                    //
                    Vector3d dir = indir.GetFaceVar(i1);
                    Vector3d bVIn = Vector3d.DotProduct(cV, dir) / dir.LengthSquared() * dir;
                    if (Vector3d.DotProduct(bVIn, face.Normal) <= 0)
                    {
                        bouU[i1] = bVIn.X;
                        bouV[i1] = bVIn.Y;
                        bouW[i1] = bVIn.Z;
                        continue;
                    }
                    //
                    Vector3d bVOut = cV + factors[i1 - start];
                    if (Vector3d.DotProduct(bVOut, face.Normal) >= 0)
                    {
                        bouU[i1] = bVOut.X;
                        bouV[i1] = bVOut.Y;
                        bouW[i1] = bVOut.Z;
                        continue;
                    }
                }
                bouU[i1] = 0;
                bouV[i1] = 0;
                bouW[i1] = 0;
            }
        }
        #endregion
        #region 更新矩阵Ab等
        public void FormAbOnFaces(ThreadInfo threadInfo)
        {
            var mbs = InSimpleFields.SimpleBounds;
            //2.调用U V W 的构造矩阵函数
            U.FormABOnFaces(threadInfo);
            V.FormABOnFaces(threadInfo);
            W.FormABOnFaces(threadInfo);
            //3.统一增加源项
            //由于粘性导致的Source
            {
                addFaceViscositySource_inner(threadInfo);
                for (int i1 = 0; i1 < mbs.Length; i1++)
                {
                    addFaceViscositySource_bound(Zone.Bounds[i1], threadInfo);
                }
            }
            for (int i1 = 0; i1 < mbs.Length; i1++)
            {
                switch (mbs[i1])
                {
                    case SIMPLEFields.NoSlipWall noslip:
                        addUVW_Ab_WallShearStress(Zone.Bounds[i1], threadInfo);//不采用fixValue的形式，可以实现层流与湍流的通用性
                        addUVW_Ab_NormalStress(Zone.Bounds[i1], threadInfo);
                        break;
                    case SIMPLEFields.Symmetry sym:
                        addUVW_Ab_NormalStress(Zone.Bounds[i1], threadInfo);
                        break;
                    case SIMPLEFields.PressureInOut p:
                        formBoundFactorOfPressreInOut(Zone.Bounds[i1], p, threadInfo);
                        break;
                }
            }
        }
        public void FormAbOnCells(ThreadInfo threadInfo)
        {
            //1
            U.FormABOnCells(threadInfo);
            V.FormABOnCells(threadInfo);
            W.FormABOnCells(threadInfo);
            //2
            {
                //附加的Source,UVW标量场中不添加,这里统一添加
                double[] U_ADiag = U.ADiag; double[] U_b = U.B;
                double[] V_ADiag = V.ADiag; double[] V_b = V.B;
                double[] W_ADiag = W.ADiag; double[] W_b = W.B;
                Vector3d[] pressureGrad = InSimpleFields.Pressure.Grad;
                (int i1, int endp1) = Zone.Inner.CellIdx.Split(threadInfo);

                if (UVW_Source_NoPressure_a.IsEqual0())
                {
                    if (UVW_Source_NoPressure_b.IsEqual0())
                    {
                        for (; i1 < endp1; i1++)
                        {
                            Vector3d sb = pressureGrad[i1] * InnCells[i1].Volumn;
                            U_b[i1] -= sb.X; V_b[i1] -= sb.Y; W_b[i1] -= sb.Z;
                        }
                    }
                    else
                    {
                        var sou_b = UVW_Source_NoPressure_b.GetCellVarArr(Zone.Inner);
                        for (; i1 < endp1; i1++)
                        {
                            Vector3d sb = (sou_b.GetCellVar(i1) - pressureGrad[i1]) * InnCells[i1].Volumn;
                            U_b[i1] += sb.X; V_b[i1] += sb.Y; W_b[i1] += sb.Z;
                        }
                    }
                }
                else
                {
                    var sou_b = UVW_Source_NoPressure_b.GetCellVarArr(Zone.Inner);
                    var sou_a = UVW_Source_NoPressure_a.GetCellVarArr(Zone.Inner);
                    for (; i1 < endp1; i1++)
                    {
                        double vol = InnCells[i1].Volumn;
                        var sa = sou_a.GetCellVar(i1) * vol;
                        U_ADiag[i1] -= sa.X; V_ADiag[i1] -= sa.Y; W_ADiag[i1] -= sa.Z;
                        var sb = (sou_b.GetCellVar(i1) - pressureGrad[i1]) * vol;
                        U_b[i1] += sb.X; V_b[i1] += sb.Y; W_b[i1] += sb.Z;
                    }
                }
            }
        }
        public void UpdateBoundsValue(ThreadInfo threadInfo)
        {
            var mbs = InSimpleFields.SimpleBounds;
            for (int i1 = 0; i1 < mbs.Length; i1++)
            {
                switch (mbs[i1])
                {
                    case SIMPLEFields.Symmetry sym:
                        setVelOfSlipWall(Zone.Bounds[i1], threadInfo);
                        break;
                    case SIMPLEFields.PressureInOut p:
                        updateBoundValueOfPressureInOut(Zone.Bounds[i1], p, threadInfo);
                        break;
                    case SIMPLEFields.NoSlipWall noslip:
                        setVelOfNoSlipWall(noslip.BoundVelocity, Zone.Bounds[i1], threadInfo);
                        break;
                    case SIMPLEFields.VelocityInlet velInlet:
                        setVelDirect(velInlet.Velocity, Zone.Bounds[i1], threadInfo);
                        break;
                }
            }
        }
        protected void addFaceViscositySource_inner(ThreadInfo threadInfo)//(Grad(Velocity))T 带来的源项
        {
            Vector3d[] U_grad = U.Grad; double[] U_bFaceCommon = U.InnFaceShare;
            Vector3d[] V_grad = V.Grad; double[] V_bFaceCommon = V.InnFaceShare;
            Vector3d[] W_grad = W.Grad; double[] W_bFaceCommon = W.InnFaceShare;
            (int i1, int endp1) = Zone.Inner.FaceIdx.Split(threadInfo);
            var gf = UVW_Gamma.GetFaceVarArr(Zone.Inner);
            for (; i1 < endp1; i1++)
            {
                ref FvInnFace face = ref InnFaces[i1];
                Vector3d dU = face.CalAverage(U_grad);
                Vector3d dV = face.CalAverage(V_grad);
                Vector3d dW = face.CalAverage(W_grad);
                double gamma = gf.GetFaceVar(i1);
                U_bFaceCommon[i1] += gamma * Vector3d.DotProduct(new Vector3d(dU.X, dV.X, dW.X), face.Normal);
                V_bFaceCommon[i1] += gamma * Vector3d.DotProduct(new Vector3d(dU.Y, dV.Y, dW.Y), face.Normal);
                W_bFaceCommon[i1] += gamma * Vector3d.DotProduct(new Vector3d(dU.Z, dV.Z, dW.Z), face.Normal);
            }
        }
        protected void addFaceViscositySource_bound(FvBound bound, ThreadInfo threadInfo)//(Grad(Velocity))T 带来的源项，理论上所有Bound都需要加上，待研究
        {
            Vector3d[] U_grad = U.Grad; double[] U_bFaceCommon = U.BouFaceShare;
            Vector3d[] V_grad = V.Grad; double[] V_bFaceCommon = V.BouFaceShare;
            Vector3d[] W_grad = W.Grad; double[] W_bFaceCommon = W.BouFaceShare;
            (int i1, int endp1) = bound.FaceIdx.Split(threadInfo);
            var gf = UVW_Gamma.GetFaceVarArr(bound);
            for (; i1 < endp1; i1++)
            {
                ref var face = ref BouFaces[i1];
                Vector3d dU = U_grad[face.OSideCell];
                Vector3d dV = V_grad[face.OSideCell];
                Vector3d dW = W_grad[face.OSideCell];
                double gamma = gf.GetFaceVar(i1);
                U_bFaceCommon[i1] += gamma * Vector3d.DotProduct(new Vector3d(dU.X, dV.X, dW.X), face.Normal);
                V_bFaceCommon[i1] += gamma * Vector3d.DotProduct(new Vector3d(dU.Y, dV.Y, dW.Y), face.Normal);
                W_bFaceCommon[i1] += gamma * Vector3d.DotProduct(new Vector3d(dU.Z, dV.Z, dW.Z), face.Normal);
            }
        }
        protected void addUVW_Ab_WallShearStress(FvBound bound, ThreadInfo threadInfo)
        {
            double[] U_ABouH = U.ABouH; double[] U_bFaceShare = U.BouFaceShare;
            double[] V_ABouH = V.ABouH; double[] V_bFaceShare = V.BouFaceShare;
            double[] W_ABouH = W.ABouH; double[] W_bFaceShare = W.BouFaceShare;
            (int i1, int endp1) = bound.FaceIdx.Split(threadInfo);
            var gf = UVW_Gamma.GetFaceVarArr(bound);
            var innVels = GetInnCellVelVars();
            var bouVels = GetBouCellVelVars();
            for (; i1 < endp1; i1++)
            {
                ref var face = ref BouFaces[i1];
                Vector3d ov = innVels.GetCellVar(face.OSideCell);
                Vector3d nv = bouVels.GetCellVar(i1);
                double XX = face.Normal.X * face.Normal.X;
                double YY = face.Normal.Y * face.Normal.Y;
                double ZZ = face.Normal.Z * face.Normal.Z;
                double XY = face.Normal.X * face.Normal.Y;
                double YZ = face.Normal.Y * face.Normal.Z;
                double XZ = face.Normal.X * face.Normal.Z;
                double commonterm = gf.GetFaceVar(i1) / Vector3d.DotProduct(face.OToN, face.Normal);
                double temp = YY + ZZ;
                U_ABouH[i1] -= commonterm * temp;
                U_bFaceShare[i1] += commonterm * (
                     nv.X * temp
                    + (ov.Y - nv.Y) * XY
                    + (ov.Z - nv.Z) * XZ);
                temp = XX + ZZ;
                V_ABouH[i1] -= commonterm * temp;
                V_bFaceShare[i1] += commonterm * (
                      (ov.X - nv.X) * XY
                    + nv.Y * temp
                    + (ov.Z - nv.Z) * YZ);
                temp = XX + YY;
                W_ABouH[i1] -= commonterm * temp;
                W_bFaceShare[i1] += commonterm * (
                      (ov.X - nv.X) * XZ
                    + (ov.Y - nv.Y) * YZ
                    + nv.Z * temp);
            }
        }
        protected void addUVW_Ab_NormalStress(FvBound bound, ThreadInfo threadInfo)
        {
            double[] U_ABouH = U.ABouH; double[] U_bFaceShare = U.BouFaceShare;
            double[] V_ABouH = V.ABouH; double[] V_bFaceShare = V.BouFaceShare;
            double[] W_ABouH = W.ABouH; double[] W_bFaceShare = W.BouFaceShare;
            (int i1, int endp1) = bound.FaceIdx.Split(threadInfo);
            var gf = UVW_Gamma.GetFaceVarArr(bound);
            var innVels = GetInnCellVelVars();
            for (; i1 < endp1; i1++)
            {
                ref var face = ref BouFaces[i1];
                Vector3d ov = innVels.GetCellVar(face.OSideCell);
                double common = 2 * gf.GetFaceVar(i1) / Vector3d.DotProduct(face.OToN, face.Normal);
                U_ABouH[i1] -= common * (face.Normal.X * face.Normal.X);
                U_bFaceShare[i1] -= common * (ov.Y * face.Normal.Y + ov.Z * face.Normal.Z) * face.Normal.X;
                V_ABouH[i1] -= common * (face.Normal.Y * face.Normal.Y);
                V_bFaceShare[i1] -= common * (ov.X * face.Normal.X + ov.Z * face.Normal.Z) * face.Normal.Y;
                W_ABouH[i1] -= common * (face.Normal.Z * face.Normal.Z);
                W_bFaceShare[i1] -= common * (ov.X * face.Normal.X + ov.Y * face.Normal.Y) * face.Normal.Z;
            }
        }
        #endregion
        public void SetZone_TransTerm(bool isEnable, IField<double> trans, IField<double> density)
        {
            U.SetZone_TransTerm(isEnable, trans, density);
            V.SetZone_TransTerm(isEnable, trans, density);
            W.SetZone_TransTerm(isEnable, trans, density);
        }
        public override void FromJsonKV(string key, IJsonReader jsReader)
        {
            switch (key)
            {
                case nameof(U): jsReader.ReadObject(U, e => U = e, t => new NSEquation()); break;
                case nameof(V): jsReader.ReadObject(V, e => V = e, t => new NSEquation()); break;
                case nameof(W): jsReader.ReadObject(W, e => W = e, t => new NSEquation()); IsVelInit = true; break;
                default: base.FromJsonKV(key, jsReader); break;
            }
        }
        public override void ToJsonKV(IJsonWriter jsWriter)
        {
            base.ToJsonKV(jsWriter);
            jsWriter.WriteObject(nameof(U), U);
            jsWriter.WriteObject(nameof(V), V);
            jsWriter.WriteObject(nameof(W), W);
        }
    }
}
//protected void updateBoundValueOfPressureInOut(FvBound bound, SIMPLEFields.PressureInOut pbound, ThreadInfo threadInfo)
//{
//    (int i1, int endp1) = bound.FaceIdx.Split(threadInfo);
//    var df = pbound.Direction.GetFaceVarArr(bound);
//    var innVels = GetInnCellVelVars();
//    var InnDc = InSimpleFields.Pressure.InnDc;
//    var innP = InSimpleFields.Pressure.InnX;
//    var bouP = InSimpleFields.Pressure.BouX;
//    var pGrad = InSimpleFields.Pressure.Grad;
//    double[] bouU = U.BouX; double[] bouV = V.BouX; double[] bouW = W.BouX;
//    double[] innU = U.InnX; double[] innV = V.InnX; double[] innW = W.InnX;
//    double[] BouMassFlux = InSimpleFields.MassFlux.BouMassFlux;
//    var factors = pbound.Factors;
//    int start = bound.FaceIdx.Start;
//    for (; i1 < endp1; i1++)
//    {
//        ref var face = ref BouFaces[i1];
//        //要区分流出还是流入
//        double massflux = BouMassFlux[i1];
//        if (massflux > 0)////由明确的流出时，确定为流出边界
//        {
//            ref Vector3d factor = ref factors[i1 - start];
//            bouU[i1] = innU[face.OSideCell] + factor.X;
//            bouV[i1] = innV[face.OSideCell] + factor.Y;
//            bouW[i1] = innW[face.OSideCell] + factor.Z;
//        }
//        else
//        {
//            //由压力驱动，在红皮书上讲的是由连续性推算，但不知道如何实现，下面可以看作是由连续性反映到压力再反映到速度上
//            Vector3d vb = innVels.GetCellVar(face.OSideCell) - Matrix3dMath.Product(InnDc[face.OSideCell], face.CalCFGradCorrectionOfBound(pGrad, innP, bouP[i1]));
//            Vector3d dir = df.GetFaceVar(i1);
//            Vector3d vbstar = Vector3d.DotProduct(vb, dir) / dir.LengthSquared() * dir;
//            bouU[i1] = vbstar.X;
//            bouV[i1] = vbstar.Y;
//            bouW[i1] = vbstar.Z;
//        }
//    }
//}
